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A novel kind of crystal order in high-density nanoconfined bilayer ice is proposed from molecular 
dynamics and density-functional theory simulations. A first-order transition is observed between a 
low-temperature proton-ordered solid and a high-temperature proton-disordered solid. The latter 
is shown to possess crystalline order for the oxygen positions, arranged on a close-packed triangular 
lattice with AA stacking. Uniquely amongst the ice phases, the triangular bilayer is characterized by 
two levels of disorder (for the bonding network and for the protons) which results in a configurational 
entropy twice that of bulk ice. 


The exceptional polymorphism of water in its crystal 
state (with 16 known bulk phases [1-4]) arises from two 
properties of its molecular bonding: Firstly, the tendency 
to form an open and flexible network, which allows for 
many different bonding topologies. Secondly, the pos¬ 
sibility of proton disorder, which results in a distinc¬ 
tion between a high-temperature disordered phase and 
a low-temperature ordered phase on the same bonding 
network. Such a transition has been observed for al¬ 
most all phases (Ih-XI, III-IX, V-XIII [2], VI-XV [3], 
VII-VIII, XII-XIV [2]). The proton disorder introduces 
a configurational entropy following the Bernal-Fowler ice 
rules [5]. Pauling’s simple estimation [6] of the entropy, 
W = (3/2)^ (where W is the number of microstates 
for N water molecules), is surprisingly accurate; indeed, 
recent numerical calculations [7] show only a very small 
variation between all proton-disordered bulk phases from 
W - 1.504^ (XII) to VF - 1.524^ (VI). 

Nanoconfined 2D water, of interest in fieids ranging 
from biochemistry to nanotechnology, has shown a simi¬ 
lar capacity for polymorphism (see review and references 
in Zhao et al [8] and Corsetti et al [9]). When con¬ 
sidered, the configurational entropy of nanoconfined ices 
has been found to be either non-extensive, ieading to an 
intrinsically ordered phase (e.g., monolayer square ice [9- 
11]), or less than or equal to that of the disordered bulk 
phases (e.g., monolayer [9] and bilayer [12] honeycomb 
ice). No order-disorder transition on the same network 
has yet been identified. 

In this paper, we investigate the phases of bilayer ice 
obtained by confining water in one dimension between in¬ 
finite plates with sub-nm separation. Different properties 
of the liquid and solid bilayer have been explored in previ¬ 
ous experimental [13] and computational studies [10-26], 
all of which have made use of empirical force-field models 
of water (mW [22], SPC/E [13, 19, 21, 25, 26], TIP4P [11, 
12, 14, 15], TIP4P/Ice [22], TIP5P [10, 16-18, 20, 23, 24], 
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FIG. 1. Schematic diagram of the stable crystalline phases 
of bilayer ice obtained in this study for a confinement width 
of 8 A. The diagram is semi-quantitative, combining enthalpy 
calculations by DFT-AIRSS with the high-temperature phase 
transition to triangular ice obtained by MD (pressures along 
the heating path are given by the thick dashed black line). 
See Ref. [27] for a more detailed phase diagram of the high- 
temperature range from MD, including the melting line which 
is omitted here for simplicity. 

ST2 [17]). The majority of studies report a low-density 
honeycomb bilayer crystal [11, 12, 15, 17, 19, 21-24, 26]. 
Han et al [23] distinguish a separate high-density solid 
(also seen by Zangi and Mark [16]), which they define as 
rhombic and suggest could be a hexatic phase. Bai and 
Zeng [24] describe a similar solid, but define it instead as 
very-high-density bilayer amorphous ice. 

We look in detail at this high-density phase, showing 
that the oxygens form a fixed and ordered triangular lat¬ 
tice which only becomes apparent when averaging over 
a sufficiently long timescale. The high degree of pro¬ 
ton disorder allowed by the lattice gives the bilayer tri¬ 
angular phase several interesting and unique properties; 
in particular, it possesses twice the entropy of the bulk 
phases without breaking the Bernal-Fowler ice rules. We 
also demonstrate the transition to at least one ordered 
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FIG. 2. The structure of triangular bilayer ice from MD at 
310 K. (a) Instantaneous snapshot, shown parallel (top panel) 
and perpendicular (bottom panel) to the confinement direc¬ 
tion. The cell side is 34.40 A. The oxygens in the two layers 
are colored differently to help visualization, (b) Oxygen po¬ 
sitions averaged over the entire run; both layers are included. 
The inset shows the average nearest-neighbor positions with 
respect to a central molecule; in this case both oxygens (in 
blue) and protons (in red) are shown. The cell side for the 
inset is 7 A. 
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FIG. 3. Potential energy per molecule as a function of tem¬ 
perature from MD. The cell shape and size and the number 
of molecules is kept constant. The insets show the structure 
of square tubes bilayer ice at 235 K and triangular bilayer ice 
at 310 K. The position of the oxygens (in blue) and protons 
(in red) is averaged over the entire run for both phases, with 
both layers being included. A (15 x 15) A^ portion of the 
entire cell is shown. 



phase at low temperature, with a possible distinct or¬ 
dered phase appearing at very high lateral pressure, as 
shown in Fig. 1. 

We use two models at different levels of theory for our 
calculations. Firstly, we perform classical molecular dy¬ 
namics (MD) simulations of 294 water molecules with the 
TIP4P/2005 [28] force-field model in the LAMMPS [29] 
code. We use a time step of 0.5 fs, a cutoff of 12 A 
for Lennard-Jones (L-J) interactions, and the particle- 
particle particle-mesh [30] (PPPM) method for long- 
range Coulombic interactions. Secondly, we perform an 
ab initio random structure search [31] (AIRSS) procedure 
to identify the lowest-enthalpy configurations for small 
unit cells (we test cells of four and eight molecules). The 
AIRSS uses density-functional theory (DFT) calculations 
in the SIESTA [32, 33] code with a fully non-local ex¬ 
change and correlation functional that gives an accurate 
description of van der Waals interactions in water from 
first principles [34]. A detailed description of the DFT 
and AIRSS methodologies can be found in our previous 
study of monolayer ice [9]. For both the MD and DFT 
simulations the confining walls are described by the same 
classical L-J 9-3 potential [9], and the simulation cell is 
periodic in all directions with a buffer region in 2 ; of 11 A. 

Fig. 2 shows a representative example of the triangular 
phase obtained by MD. The simulation is carried out in 
a fixed square cell with an area per molecule of 4.02 A^ 
and a confinement width of 8 A. The cell is first equili¬ 
brated with a Berendsen thermostat for 60 ns, after which 
statistics are collected for 100 ps in the NVE ensemble. 

Any instantaneous snapshot of the system (Fig. 2(a)) 
seems to confirm the definition of an amorphous solid. 
Although the particles are densely packed, they are not 


arranged on a clear lattice and there is no visible order 
in the hydrogen-bonding network; instead, there is an 
irregular combination of squares, pentagons, hexagons, 
and so on (similarly to low-density amorphous bilayer 
ice [11, 15, 17, 24]). If quenched to low temperature 
the system would therefore freeze into this amorphous 
arrangement. However, the bonding network constantly 
rearranges during the course of the simulation, and the 
molecules are pulled towards different neighbors depend¬ 
ing on the instantaneous bonding arrangement. When 
averaging the oxygen positions over the entire 100 ps 
(Fig. 2(b)) a very clear triangular lattice emerges. The 
unit cell is found to be hexagonal to a good approxima¬ 
tion (with errors of 1% in the unit cell angle and 0.3% 
in the a/h ratio, which can be attributed to finite size 
effects due to the fixed supercell). The inset of the figure 
shows that every molecule is at the centre of a hexago¬ 
nal cage with six in-plane neighbors; the hydrogen bonds 
(shown by the inner circle of protons bonding to the cen¬ 
tral molecule) fluctuate between all six neighbors with no 
preferential direction. It is also important to note that 
the bilayer has an A A stacking. Each molecule therefore 
maintains a constant bond to the other layer, and three 
fluctuating in-plane bonds. The number of molecules 
with defective configurations is reasonably small, <5% 
of the total. 

The molecular bonding arrangement of bilayer trian¬ 
gular ice is consistent with the expected sixfold orienta¬ 
tional order of a hexatic phase; we discuss this in more 
detail elsewhere [27]. It is interesting to note that a 
hexatic ice monolayer has previously been predicted and 
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Bond configurations 



FIG. 4. All non-equivalent hydrogen-bond configurations 
around a site of the triangular lattice (symmetry degeneracy 
given in brackets). Some examples of extended networks are 
shown below; for the square tubes and square wave networks 
the bonds in the second layer are shown in light blue. 


analyzed using a coarse-grained model [35]; however, in 
that case the number of hydrogen bonds was found to be 
greatly reduced, and, therefore, is significantly different 
to what we observe. 

The fluctuating hydrogen bonds suggest a large config¬ 
urational entropy and the possibility of an order-disorder 
transition upon cooling. We first verify this empiri¬ 
cally within our MD simulation, by starting from a high- 
temperature configuration (prepared as described previ¬ 
ously) and cooling in steps of 5 K. Each step is carried 
out by re-equilibrating the system with a Berendsen ther¬ 
mostat for 5 ns, and then collecting NVE statistics for 
100 ps. The reverse heating process is also performed 
independently. The two processes are shown in Fig. 3. 

The sharp discontinuity and hysteresis loop in the po¬ 
tential energy of the system are indications of a first-order 
phase transition occurring at Tc — 280 ± 15 K. Below 
this temperature the triangular bilayer transforms into a 
phase characterized by rows of square-based tubes, with 
no hydrogen bonds between them [24]. The unit cell 
symmetry is lowered from hexagonal to rhombic (cen¬ 
tred rectangular), with an angle of ^106°. The posi¬ 
tion of the protons as well as the oxygens in the square 
tubes is fixed for the entire simulation, denoting a proton- 
ordered phase with no configurational entropy. We can 
therefore estimate the configurational entropy of the tri¬ 
angular phase by equating it with the gain in inter¬ 
nal energy at the transition (assuming a similar vibra¬ 
tional contribution). Since we are in the NVE ensemble, 
NS = NU/Tc ^ 0.7 ± 0.1 /cb/ molecule. The configu¬ 
rational entropy is therefore almost twice the value for 
bulk ice, S/N 0.4 /cb- Additional calculations in the 
NPxyT ensemble give a perfect qualitative and quantita¬ 
tive agreement both for Tc and AS (estimated in this case 
as AH/Tc); details are in the Supplemental Material. 

The large entropy value can be understood by consid¬ 
ering the generalization of the ice rules to the case of the 
triangular bilayer lattice. There are two levels of config¬ 
urational entropy, in the arrangement of hydrogen bonds 
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a molecule on the triangular bilayer lattice (symmetry de¬ 
generacy given in brackets). Out-of-plane molecules donate a 
proton to the interlayer bond (not pictured). The frequency 
for each configuration, shown below, is calculated from MD 
simulations of the triangular phase. 


on the lattice and in the placement of protons along those 
bonds. For all other solid phases only the latter applies, 
since the bonding network is fixed. 

Fig. 4 shows the possible configurations of the three 
in-plane hydrogen bonds around a molecule. As already 
mentioned, the fourth bond is fixed and connects the 
molecule to the other layer. The bonding pattern in the 
two layers is therefore decoupled. From the random net¬ 
work example shown in the figure it can be seen how 
polygons of varying size and shape are created, result¬ 
ing in the characteristic amorphous appearance of the 
instantaneous snapshots. 

It should also be expected that low-energy configu¬ 
rations will result from particular regular bonding pat¬ 
terns; indeed, the square tubes network can be created on 
the triangular lattice entirely from B configurations (see 
Fig. 4). This is an unusual pattern, made up of entirely 
disconnected networks with ID periodicity. The regular 
pattern allows for a macroscopic distortion of the lattice 
(as opposed to the localized distortions of the molecules 
in the random network), which is enhanced by the weak 
interaction between tubes. This results in the change of 
symmetry observed in the simulation. 

The second level of configurational entropy is the ar¬ 
rangement of protons on the bonding network (Fig. 5), 
such that each molecule has two short and two long OH 
distances. The two layers are now no longer indepen¬ 
dent, because the total number of in-plane and out-of¬ 
plane molecules in the system has to be equal, but they 
do not have to be equally distributed between the two 
layers. There are 12 non-symmetrically equivalent pro¬ 
ton configurations for an individual molecule, and 120 
distinguishable configurations in total. 
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FIG. 6. Relaxed crystal structures for stable bilayer ice phases 
calculated by DFT and found with an AIRSS methodology, 
(a) Honeycomb, lateral pressure of 0.01 GPa-A, area per 
molecule of 4.98 A^. (b) Square tubes, lateral pressure of 

0.01 GPa-A, area per molecule of 4.03 A^. (c) Square wave, 
lateral pressure of 10 GPa-A, area per molecule of 3.83 A^. 
The confinement width is 8 A. The four-molecule unit cell is 
shown. 


It should be expected that not all proton configura¬ 
tions are energetically equivalent, and some might be 
effectively prohibited for this reason, thus reducing S. 
Table SI in the Supplemental Material gives our Pauling- 
like estimate for the configurational entropy of the sys¬ 
tem depending on which proton configurations are al¬ 
lowed; we use a Monte Carlo method to calculate this, 
also detailed in the Supplemental Material. Including all 
possible configurations gives an unreasonably high value 
of 1.35 ± 0.01 /cb- Indeed, if we analyze our MD sim¬ 
ulations, as shown in the histogram in Fig. 5, we find 
that four configurations are almost entirely absent (53, 
Cl, C2, C2*). Interestingly, there is also a noticeable 
asymmetry between in-plane and out-of-plane configura¬ 
tions, with the latter being more permissive (53* and 
Cl* are allowed). Taking these restrictions into account, 
the Pauling estimate is 0.8 ± 0.1 /cb, in good agreement 
with the value obtained previously from MD. 

We conclude our analysis by returning to the overall 
phase diagram presented in Fig. 1. Here we make use 
of first principles calculations with DFT and the AIRSS 
methodology to identify the lowest-enthalpy structures 
at different lateral pressures. Temperature dependence 
is obtained by estimating the phases’ configurational en¬ 
tropy, and, hence, Gibbs free energy. Our previous study 
of monolayer ice has shown the contribution to free en¬ 
ergy differences from vibrational effects to be small [9]. 

The AIRSS recovers the low-density honeycomb bi¬ 
layer crystal (Fig. 6(a)). Contrary to previous re¬ 
ports [12], there is no noticeable distortion in the hexag¬ 
onal unit cell, and so no preferential proton orientation. 
The derivation of the configurational entropy is there¬ 
fore equivalent to that of Pauling for bulk ice (more than 
twice the previous estimate for the bilayer [12]). 

Two high-density bilayer crystals, close in enthalpy, 
are predicted: square tubes ice (Fig. 6(b)), and a novel 


structure not found by MD which we refer to as square 
wave ice (Fig. 6(c)). The latter has a similar bonding 
network to the tubes, but with the two layers shifted rel¬ 
ative to each other; the whole crystal is now connected 
into a single network, and when viewed perpendicular to 
the confinement direction gives a square wave pattern. 
This is illustrated schematically in Fig. 4. The rhombic 
distortion with respect to the hexagonal cell is less pro¬ 
nounced for the wave than for the tubes, with unit cell 
angles of 112° and 101°, respectively. 

The square wave and the square tubes phases are pre¬ 
dicted to be ordered on both levels. This is because the 
bonding network is fixed, and all the examples recovered 
from the AIRSS make use of only four proton configu¬ 
rations (51, 51*, 52, 52*). Two other possible con¬ 
figurations compatible with the network are not found 
(53, 53*). It can easily be shown that this results in a 
non-extensive entropy. It is important to note that the 
specific configuration of protons in the square tubes ob¬ 
tained by MD using TIP4P/2005 is in agreement with 
the lowest-energy configuration found by DFT. A previ¬ 
ous study using TIP5P instead gives a different configu¬ 
ration [24]. This is a similar result to what was found for 
the square monolayer [9], and a further indication of the 
accuracy of the TIP4P/2005 force-field model. 

Recent experimental observations of ice confined be¬ 
tween graphene sheets [13] show a square configuration 
for the monolayer up to the trilayer. While DFT and 
some force-field models are in excellent agreement with 
experiment for the monolayer [9], we have been unsuc¬ 
cessful in obtaining the square bilayer; calculations using 
both DFT and TIP4P/2005 find it to be unstable with 
respect to the high-density phases discussed here. This 
suggests interesting effects originating either from deeper 
levels of theory, or considerations such as the confinement 
and the dynamics of formation, meriting further study 
both from experiment and simulation [36]. 
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